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Abstract 


This paper presents a computational method for solving two types of 
integro-differential equations, system of nonlinear high order Volterra-Fredholm 
integro-differential equation(VFIDEs) and nonlinear fractional order integro- 
differential equations. Our tools for this aims is operational matrices of inte- 
gration and fractional integration. By this method the given problems reduce 
to solve a system of algebraic equations. Illustrative examples are included 
to demonstrate the efficiency and high accuracy of the method. 


Keywords: Operational matrix of integration; Volterra-Fredholm; Non- 
linear system of integro-differential equations; Fractional order; Legendre 
wavelet. 


1 Introduction 


Integro-differential equations frequently appear in all fields of sciences such 
as physics, chemistry and engineering problems [11, 20, 23, 24]. In last few 
decades fractional calculus and fractional differential equations have found 
application in several different disciplines, many important phenomena in 
electromagnetic, acoustics, viscolasticity, electrochemistry and material sci- 
ence are well described by differentiable and integro differentiable equation 
of fractional order[3, 22]. There are various numerical and analytical meth- 
ods to solve such problems, for example, the homotopy perturbation method 
[4, 7, 8, 9], the Adomian decomposition method [5], fractional differential 
transform method [21] and Gronwald-Letnikov discretization method [6]. 
In recent years the approximation of orthogonal functions has been play- 
ing role in the solution of different kinds of mathematical and engineering 
problems such as identification, analysis and optimal control[{15, 16, 18]. The 
main feature of this technique is to reduce the integro-differential equations 
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to a nonlinear algebraic equation by introducing integration matrix of basis 
functions. In present article, we are concerned with the application of Leg- 
endre wavelet to the numerical solution of: 

(I). Nonlinear fractional order integro-differential equations 


D?,u(t) a+ [ K(t, u(t), Do, u(t)) dt, O<a<l. (1) 
(II). Nonlinear system of high order (VFIDe) of the form 


> pij(x)u (x) = fila) + a fe’ Ka(2,t, u(t), w/(t),...,u™ (tat 
+Ni2 [i K io(x, t, u(t), u/(t),...,u™(#))dt, i =1,...,8, 


where uJ) (x) = Cae ae ug?) for 7 = 0,...,n and initial conditions are 


u(0) =a;, 7 =0,1...,n—-1, (3) 


where f(x),K, Kj, and Kjg are known functions assumed to be in L?(R) 
on the interval 0 < x,t < 1, u(t) is unknown, Kj, and Kjz are nonlinear in 
a,t,u(t),...,u((t). This type of equations whose integrand contain high 
order derivatives arise in many fields such as theory of elasticity . 


The article is organized as follows: in Section 2 we define the Legendre 
wavelets and operational matrix of integration. Section 3 is devoted to the 
solution of Eq. (1). In Section 4, we obtain an error bound for our method. 
Section 5, include our numerical findings and demonstrate the accuracy of 
the proposed scheme. 


2 Preliminaries and notation 


This section gives some necessary definition and mathematical preliminaries 
of the fractional calculus theory which are used further in this paper. The 
Riemann-Lioville fractional integration of order a > 0 is defined as [14] 
ii 
tee (¢) ~ Ta) i G— oe a) dr, (4) 
P£(t) = f(r), 


and its fractional derivative of order a > 0 is normally used: 


d 


Dif) = (Ge F() (n-1<asn), (5) 


where n is an integer. For Riemann-Lioville definition, one has 
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T(u + 1) 
Tee? = te, 6 
: T(a+u+1) (6) 


The modified fractional differential operator D%, proposed by Caputo is 
t n—-a-1 p(n 
ranay lo (t—T) f' )(r) (n-1l<a<n), 


pas=4 (7 
“ f(t) a=neNn, 


where n is an integer. Caputos integral operator has an useful property: 


n-1 & 


Ie D2, f(t) = f() - YOON (n-1<a<n), (8) 


where n is an integer. 


3 Properties of Legendre wavelet 


3.1 Wavelets and Legendre wavelet 


Wavelet constitute a family of functions constructed by a single function 
called the mother wavelet. When the dilation parameter a and translation 
parameter b vary continuously, we have the following family of continuous 
wavelet as [10] 


= i=5 
dia) = lal Ue) abe R,aF0. 


If we restrict the parameters a and b to discrete values as a = aj’”,b = 
kbpag”, ao > 1,60 > 0 and m,k € Z.We have the following family of discrete 
wavelets 
Vm.k(t) = |aol""/*b(ag't — kbo), 
where Wm.z(t) forms a wavelet basis for L?(R). In particular, when ag = 2, 
bo = 1, Wm.«(t) forms an orthonormal basis. 
The Legendre wavelets are defined on interval [0,1) see [16, 17]. 


4/m+ 52/2 Din (2*t — fn), for=4 <t< am 


0 otherwise, 


Vnm= 


where m = 0,1,...M—1 and n=1,2,3,...,2*~!. The coefficient ,/m + 3 
is for orthogonality. Here, L,,(t) are the well-known Legendre polynomials 
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of order m which are defined on the interval [- 1,1] and can be determined 
with the aid of the following recurrence formulae: 


Lo(t)=1, Lilt) =t, 


2m+1 m 
thm (t) — (——~) Lm_ : S128 0.s 
tml t) = (> )Em-r(t), 


3.2 Function approximation 


Theorem. A function f(t) defined on [0,1) can be expanded as infinite sum 
of Legendre wavelets, and the series converges uniformly to the function f(x), 


that is a act 

0 
where, Cnm = (f(t), nm(t)), in which (.,.) denote the inner product. 
Proof. see[13]. 


If the infinite series in Eq. (9) is truncated, then it can be written as 


2k-l M1 


[D2S. S Gant =O Ue), (10) 


n=1 m=1 


where C and W(t) are 2” x 1 matrices given by 


T 
C = (C10, C11, +--+, CLM—1; C20,+- +, C2M—1,+++5 Cok-19,--+ 5 Cok-1yy] , (11) 


W(t) = [v10, Wu, ere »WimM-1) P20; eae ,WeM—1; eee Wor-19, See ,Wor-im)’. 
(12) 
Now we want to find an upper bound to the estimate error . Suppose that 
f(x) is a (m+ 1)—times differentiable function on Q = [0,1). An error 
function between f(x) and its Legendre-wavelet approximation fnm(x) is 
defined on every subinterval Q,, = [434 <t < 4] as 


€nm(X) = f(x) — fam(x) = f() — CrnmYnm(z). (13) 


Then we can write 


ll enum (a) =f |f(2) — Cam nm(a))’. (14) 


ok 


Since Unm(x) is a polynomial of degree m,we can use the error bound for 
interpolation of degree m on Q,, that is 
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pim+)) 
If(@) — pala) < A max |g], (15) 

4(m + 1) ec[a2 a] 

where h = 5% Eq. (14)and Eq. (15) 
2 
2 oe | amy (m+1) 
Venm (a) < faze fam 5 Ma HMO 
a 2: (16) 
(m+1) mn 
S ol dqmaty ,_ max, [FOP 


ge[ 252, 43] 


According to above equation we find an error bound for each subinterval as 


1 Alm) was i 
lnm (I< 5579 aa T) cep (6)|. (17) 
Then for error on 2 we get 
1 A(mt1) 
eta forty) 
lel < ee ameD Bey OPO}: (18) 


3.3. The Legendre wavelets operational matrix of 
integration 


The integration of the Vector defined in Eq.(12) can be obtained as 


eG W(t')dt’ = PW(t), (19) 
0 


where P is the 2'-!M x 2*-!M operational matrix for integration [18] 


LHHH.-H 
OLHH-:--H 
,|OO0OLH.-.H 


o°o 


Org ee i TE 
00 0::. L 


A and L are M x M matrices given by : 
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2 0.:-- 0 
0 0O.:-- 0 
H= ; 
0 0 .:-- 0 
and 
| 1 ats 0 0 0 0 0 
3i/? 3i/? 
ae} 0 3x 51/2 0 0 0 
51/2 51/2 
0 saa 0 ene 0 0 0 
_ 7/2 
L= 0 0 xez «(COO 0 0 0 
. : ; , (2M—3)'/2 , (2m_3)1/2 
0 0 0 Ol se ~ (QM —3)(2M—5)1/2 0 (2M—3)(2M—1)!/2 
0 0 0 0 0 See ae 0 
ae ~ QM—1)(2M-3)1/2 


3.4 Operational matrix of fractional integration 


We defined a m-set of Block Pulse function (BPF)as: 


Hy dy LES EEL) pt 
b(t) = 5 otherwise, (20) 
where 7 = 0,1,2,...(m— 1). 
The function 6;(t) are disjoint and orthogonal. That is 
0) #43 
(noi ={ Py, Per (21) 


The Legendre wavelet may be expanded into m-terms of block pulse function 
(BPF) as 
where ;, 
Bm (t) = [bo(t) bi (t) ... bi(t) ... bin-ay@]}": (23) 
The Block Pulse operational matrix of the fractional integration give in [12] 
F° as following: 
(IP Bm)(t) ¥ F° Bm(t), (24) 


where 
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Lei Corte G ii) 

F , ODS Sa) 

PEO a ee OL E(m-—3) 
mo (a+ 2) ; : 
000 °-. : 
000 0 1 


with & = (k+1)%*!—2k°t!+(k—1)°*!. The Legendre wavelet operational 
matrix of of fractional integration is defined in [19] as 


i nee Omxml : (25) 
so the fractional integration of vector in Eq. (12) is defined as 
(pW) (t) & P°W,, (t). (26) 


4 Application to nonlinear system of VFIDEs 


Here, before presenting our method, we prove the next lemma. By this 
lemma we can approximate the high order derivative of a function by Legen- 
dre wavelet. 


Lemma. Suppose that u(x) = C?W(x) where C and V(x) are defined 
in Eq. (11) and Eq. (12), then 


u*) (2) = (CTP Ys ) BT PiU (x), (27) 
where P is operational matrix of integration, u(0) = u and E is defined 
as E7 W(t) =1 


Proof. suppose that f(z) = u*(x) and we approximate u(x) and f(x) b 
Legender wavelet as 


u(z) = CT W(z), 
{ fe) = FTY(a), i 
by integrating f(t) on (0, t] 
: f(t) dt! ...dt? f. FTP W(t!) dt’... dt! 
) i‘? a i at 0 ee 
pana times (k—1)times 
— tt t mT p2 ! ! ’ 
=foloowdf Pee) a... dt (29) 
(k—2)times 


= FT PFW(t), 
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Since f(z) = u*(zx), then 


Jo So So £8) eta dtl = fo So So (u® Yu Y(0)) aCe!) «+ at! 


k—times (k—1)—times 


= t pt t 
= fo io fo uN) (4) dt! ---dt! ul Y(0) fp Joc so dts: :dt 


(k—1)—times (k—1)—times 


= u(t) — wu (0) oad u) (0) aa all Y(0 (0) fo Soo fo did t 
veraret 1)—times 


| (30) 
by u(0) = ul we get 


uli Ae faa =u ET P(t), (31) 


k—times 
Eq. (29)-(31) result 
FT PFW(t) = CU(t) — uO ET W(t) — UY ET P(t) — --- — ul ET PE (Et) 
k-1 : 
= CTW(t) — Yi uM ETP U(E), 
i=0 
(32) 


Since the basis functions are linear independent, we omit W(t) from both 
sides of Eq. (32), then this equation can be written as 


k-1 


PPO ey iy err, (33) 

i=0 

and then re 
TOF p=k = se ul ET pi-k (34) 

according to Eq. (28) = 
u*) (2) =(CTP-* -¥ u) ET Pi’) U(a). (35) 


This ends the proof of lemma. 


To solve Eq. (2) by Legendre wavelets, we assume that each ue() has the 


expansion as 
wpa) SCP Ve), ba Ays is (36) 


by Eq. (35) the derivative expansion is given by 
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k- 
ys”) (x) = (CP P-* 2546 DET Pi-)G(x), €=1,...,8, (37) 
1=0 


substituting Eq. (36) and Eq. (37) in Eq. (2) results for @=1,...,s 


peoC} V(x) + 3 pe(x)(CPP-* - pas ul BT Pm—i)yU (a) = fo(x) 


i=l 


n=l 4. . 
tra fy Kei(a,t, CP U(t),..., (CTP — DY uM ETP") W(a)) dt 


i=0 
n—-1 . : 
+e2 fo Keo(a,t, CP U(t),..., (CTP? — Tu ET PH) W(x) at. 
i=0 
(38) 
by suitable collocation points, the zeros of Chebyshve polynomials [16] 
2i-1 
oj = cos( EU) ¢=1,...,2°-'M, (39) 


we collocate the Eq. (38). In order to use the Gaussian integration formula 
for Eq. (38), we transfer the t-intervals [0, x;] and [0,1] into ¢, and ¢2 intervals 
[-1,1] by 


2 
q= to 1 Cg = 2t— 1. (40) 
Let 


n-1 |. F 
Ha(2;,t) = Ka(a;,t, CP V(t),...,(C7P-” — XY yP ETP") U(a)), 


i=0 


Heo(a;,t) = Kea(a;,t, CP V(t),...,(C7P-” — > u ET P*-") U(x), 
i=0 
(41) 
We rewrite Eq. (38) as 
peoC7? (x5) + % pala; (CPP — = ul) ET P™-)U(a;) = felay) 


ad 2 am Ae (25,3 = (¢, + 1))ddi 
+42 ft Hala; i(Gt))de, €=1,...,8, 


(42) 
and with the Gaussian integration 


peoC? U(x5) + % pales (CF ti > ul) EP P™4)U (ars) = fela;) 


m=0 


+r 3 wines (x5, (Cin + 1)) 
h=l 


82 
geo won Heo(x;, $(Con + 1)), £= Tyne 8, 


(43) 
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where (i, and Coz are 5s; and sg zeros of Legendre polynomials L,,41 and 
L.,41 respectively, and w 1, wg; are the corresponding weights. If we assume 
that 


a—1 


Ae(e) = peoC? Was) + Y pei(es (CEP ~ ¥ ule ETP”) (x5) 


m=0 
ok Ls $2 
ro YS winHer (tj, 2(Gan + 1) — 42 YS wonHee(2j,4 (Gr + 1), 
h=1 


h=1 
Be(x) = fe(x), | te ER 2 
(44) 
Then our problem has the next matrix representation form 

Ai (21) By (21) 

Aj (&gx-1 7) By (xge-1y4) 
As(#1) B,(#1) 

Ag(Xgx-1 yy) Bs (gx-1y) 


This 2*-! Ms x2*-! Ms nonlinear system of equations which can be solved 
using Newton iterative method for the elements of C. 


5 Application to nonlinear fractional order 
integro-differential equations 


In this section we want to apply the operational matrix of fractional in- 
tegration to fractional order integro-differential equation. Assume that we 
approximate D%,u(x) by Legndre wavelet as 


Dou(z) = KTU(2), (45) 
then Eq. (8)and Eq. (26) result 
u(x) = KT P%,.,,, U(x) + u(0). (46) 


By Eq. (45)and Eq. (46) we rewrite Eq. (1) as 


KT W(x) = f(z) +f k(t, KT P&,, W(t) + u(0), kK? U(t))dt (47) 
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Assume that 


H(t) = k(t, K7 Prim Y(t) + u(0), K7 Ud), (48) 
and like the last chapter, we collocated this equation by Eq. (39)in 2*-!M 
points and then use the Gaussian integration. Finally, we can write Eq. (39) 
as 


KTW(z;) = f(2:) Se one (Cn +1)) i=1,...,2* 1M (49) 
h=1 


which is the 2*-!M x 2*-!M nonlinear of system equation which can be 
solved using Newton iterative method for the elements of C’. 


6 Numerical examples 


In this section we consider some examples which show that operational ma- 
trices are powerful and demonstrate the accuracy of our method. 


Example 5.1. Consider the nonlinear system of integro-differential equa- 
tion 


3auy (x) + uf (x) = 5a + 2uh(a) — fy (uh(t) + ur (t)uY (t))dt, + fo eae Jul, (t)dt, 

2us(x) + uf (a) = —4x? — ee r+ fo (taus(t)uy (t ) + us(t) )dt + fo a x?ug(t) + uh (t)uf (t)dt 
#/8ys(x) + 05 ( 2) =2— $a + uf? (x) — Qut(a) + fy (w?ua(t) + uw? (t) + Pug (t) jdt + fo x?ul (t)dt 
u1(0) = u4(0) = 0, u2(0) = 0, u(0) = 1, ug (0) = ug (0) = 0, 


(50) 
which has the exact solution ui(x) = x”, ue(x) = # and us(x) = 3x7. Fig- 
ure.1 show the absolute error when we apply our method for M = 3 and 
k=1. 

It is clear form figures that our approximate solution is in good agreement 
with exact one. 


Example 5.2. As a second example, consider the nonlinear system given 
in [2, 1] 


ae =1- $u(x) + fete x — t)ug(t) + uz (t)ue(t)) dé, 
ub(x) = 2x lh ( ((x — t)uy(t) — u2(t) + u?(t))dt, (51) 


which has the exact solution u;(x) = sinh(x) and u2(x) = cosh(x) for M = 6 
and k = 1. 
Results for Example 5.2 are reported in Table 1 for wi(x;) and u2(a;). 
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0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 


02 
(a). Error For u4(x) (b). Error For u(x) 
1.5107 + 
1.x 10-1 F 
5.x 10716 - 
*00 02a 10 


(c). Error For u3() 


Figure 1: Absolute error for Example 5.1 for M=3 And k=1 


Example 5.3. Consider the nonlinear fractional order integro differential 
equation given in[7| 


D%u(t) =1+ | u(t)D&u(t)dt O<a2<l0<a<1 (52) 
0 


The exact solution of this problem for a = 1 is V2Tan(2t) we solve this 
equation for m = 20 and different a numerical results are shown in Figure 2. 


Example 5.4. Finally Consider the nonlinear fractional order integro dif- 
ferential equations in[7] 


Deyu(t) = —1+ f u2(t)dt0<x<10<a<l (53) 


subject to the initial conditions y(0) = 0. Table 2 shows the numerical 
results for a = 0.8,0.9,1 when m = 20. From Table 2 we can see that 
the approximate solutions obtained by our method are in good agreement 
with the exact solution for a = 1, and with the approximate solutions for 
a = 0.8, 0.9 in [7]. 
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Table 1: Numerical result of Example 5.2 


Error for Error for Method of [2] Error for Method of [1] 
M=6,k=1 N=5 N=6 
Bi u2(x) u1(x) Uu2(x) u1(2) u2(x) u1(x) 
0  1.70x10-® 7.66 x 107 0 0 0 0 
0.1 6.82x1077 3.42 x 107 1.3 x 10-8 1x 10-8 1.41 x 107° 1.41 x 107° 


7.98x 10-7 1.331077 9.15x107-8 9.15 x 1078 
9.06 x 10-& = 2.17x107-& 1.06 107° 1.06 x 107° 
5.06 x 10-5 1.53x 10-5 6.03x 107° 6.03 x 107§ 
1.90x 1074 6641075 2.341075 2.34 1075 
5.05 x 10-* 2.12 10-4 7.08x10-5 7.08 x 10-5 
1.36 x 10-3) =5.27x 107-4 1.811075 1.81 x 10-5 
2.87x 10-8 1.05x 10-3 410x10-* 4.10 x 10-4 
5.34x 10-3 1.66 10-3 8.45x10-4 8.45 x 10-4 
8.71 x 10-3 1.17x 10-3 =1.62x 10-3 1.62 x 10-3 


0.2 285x107 8.94~x 
0.3 5.17x 10-7 2.50 x 
0.4 1117x1077 1.22x 
0.5 5.431077 2.40x 
0.6 1.751077 1.08 x 
0.7 465x10-7 2.43 x 
0.8 2.681077 4.02 x 
0.9 7.65x10-7 5.10x 

1 2.85x10-8 1.44 x 


SS SS See ee aes 
@AAIAIAN WN DAA 


Figure 2: Numerical result for Example 5.3 for different a and m=20 


7 Conclusion 


Most nonlinear integro-differential equation with nonlinear differential part 
are usually difficult to solve analytically. In many cases it is required to 
obtain the approximate solution. We have shown that the properties of op- 
erational of matrix of integration and operational matrix of fractional inte- 
gration together with Legendre wavelet can reduce the system of nonlinear 
integro-differential equation and nonlinear fractional order integro differential 
equation to a system of algebraic equations. The advantage of this method 
is that it can solve high and fractional order integro-differential equation eas- 
ier and more time efficient. Also we found an error bound. Although we 
solved our problem by Legender wavelet, other orthogonal basis also can be 
used. Illustrative examples show the high accuracy of the method in compar- 
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Table 2: Numerical result of Example 5.4 
ty Exact solution a= 1 Our method Method of [7] 
aad a=09 a=0.8 a=1 a=09 a=0.8 
0 0 0 -0.00013 — -0.00046 0 0 0 

0.0625 -0.06250 -0.06249 = -0.08574 = -0.11683 -0.06250 -0.08574 -0.11682 
0.125 -0.12498 -0.124977 -0.15995 — -0.20327 -0.12498 -0.15997 -0.20328 
0.1875 -0.18740 -0.18749 — -0.23023 — -0.28080 -0.18740 -0.23024 —-0.28082 
0.2500 -0.24968 -0.24966  -0.29788  -0.35269 -0.24968 -0.29790 -0.35272 
0.3125 -0.31171 -0.31172 — -0.36339 — -0.42026 -0.31171 -0.36342 — -0.42039 
0.3750 -0.37336 -0.37333  -0.42695 — -0.48409 -0.37336 -0.42689  -0.48413 
0.4375 -0.43446 -0.43443 — -0.48858 —-0.54446 -0.43446 -0.48861  -0.54451 
0.5000 -0.49482 -0.49478  -0.54818 — -0.60140 -0.49482 -0.54824 -0.60150 
0.5625 -0.55423 -0.55418 -0.60565 = -0.65501 -0.55423  -0.60571 -0.65510 
0.6250 -0.61243 -0.61237 — -0.66078 -0.70511 -0.61243 -0.66086 -0.70521 
0.6875 -0.66917 -0.66910 —-0.71337  -0.75162 -0.66917 -0.71345  -0.75172 
0.7500 -0.72415 -0.72418  -0.76318  -0.79440 -0.72415  -0.76327  -0.79451 
0.8125 -0.77710 -0.77710 —-0.80997 — -0.83330 -0.77710 -0.81006  -0.83341 
0.8750 -0.82767 -0.82771  -0.85348 — -0.86820 -0.82767 -0.85395 — -0.86831 
0.9375 -0.87557 -0.87564  -0.89349  -0.89896 -0.87557 -0.89361 -0.89908 


ison with other methods. This procedure can also be used for solving other 
functional equations such as ordinary and partial differential equations. 
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